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Abstract 

We derive a master equation for the electron transport through molecular wires 
in the limit of strong Coulomb repulsion. This approach is applied to two typical 
situations: First, we study transport through an open conduction channel for which 
we find that the current exhibits an ohmic-like behaviour. Second, we explore the 
transport properties of a bridged molecular wire, where the current decays exponen- 
tially as a function of the wire length. For both situations, we discuss the differences 
to the case of non- interacting electrons. 
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1 Introduction 



In recent years, it became possible to adsorb organic molecules via thiol groups 
to a metallic gold surface and, thus, to establish a stable contact between the 
molecule and the gold. This opened the way to reproducible measurements of 
the current through single molecules. Experiments for such molecular conduc- 
tance can be achieved in essentially two ways: One possible setup is an open 
break junction bridged by a molecule [1-3]. There, the current measurement 
provides evidence for single molecule conductance because asymmetries in the 
current- voltage characteristics reflect asymmetries of the molecule [3,4]. Alter- 
natively, one can use a gold substrate as a contact and grow a self-assembled 
monolayer of molecules on it. The other contact is provided by a gold cluster 
on top of a scanning tunnelling microscope tip which contacts one or a few 
molecules on the substrate [5,6]. 
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Fig. 1. Tight-binding model for an open conduction channel with N = 5 sites. An 
external bias voltage V = (/Ur. — is applied to the molecular wire. 

Naturally, the experimental effort with such molecular wires is accompanied 
by a vivid theoretical interest [7-9] . Presently, the main theoretical focus lies 
on the ab-initio computation of the orbitals relevant for the motion of excess 
charges through the molecular wire [10-14]. 

Another line of research employs rather generic models to gain a qualitative 
understanding of the transport mechanisms involved. The treatment of these 
models can be distinguished according to the level at which interaction is taken 
into account. Here, we are in particular interested in two cases in which the 
many-body problem can be traced back to the dynamics of single electrons 
on the wire: The first case premises non-interacting electrons for which the 
current can be computed from a Landauer-like formula [15-20]. The second 
case deals with the opposite limit in which Coulomb repulsion is so strong 
that at most one excess electron can be located on the molecule. Such theories 
have been developed in the context of conduction through coupled quantum 
dots [21-23] and for the incoherent transport through molecular wires [24-26]. 

In this work, we derive a master equation approach for molecular conduction 
in the limit of strong Coulomb repulsion which restricts the population of the 
molecular orbitals to zero or one excess electron. Thereby, particular care will 
be taken in avoiding inconsistencies like spurious non-vanishing transport in 
equilibrium situations. We present in Section 2 our working model and derive 
in Section 3 a master equation which we evaluate for the two mentioned limits, 
namely non-interacting electrons and strong Coulomb repulsion. Subsequently, 
we study in Sections 4 and 5 transport through open conduction channels and 
across bridged molecular wires, respectively. Explicit analytical expressions for 
a wire that consists of only two sites are derived in Appendix A. 



2 Model 



The setup at hand for studying coherent quantum transport is depicted in 
Fig. 1. The corresponding Hamiltonian reads 

H H w [ rc -\- -/^lcads H wl] - C — lead; (1) 
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where the individual terms describe the molecular wire, the electron reservoirs 
of the leads and the coupling of the wire to the leads. The wire itself is treated 
in a tight-binding approximation consisting of iV orbitals. Since we aim at 
exploring blocking effects, the corresponding wire Hamiltonian, incorporating 
the Coulomb repulsion in the limit of a large interaction strength U, assumes 
the form 

N-l 

H wirc = ]T E n clc n -A^ (c[ +1 c n + cic n+1 ) + UN(N - 1). (2) 

n n =l 

The fermion operators c\ (c n ) create (annihilate) an electron in the orbital \n) 
and E n denotes the respective on-site energy. Here, we neglect the influence 
of the voltage profile on the on-site energies [27,28]. In the Coulomb interac- 
tion term, N = J2n c li c n is the operator counting the excess electrons on the 
wire. The inter-site coupling characterised by the hopping matrix element A 
is assumed to be equal between all neighbouring sites. The leads attached to 
the molecular wire are modelled by ideal Fermi gases, 

-Pleads = € q C lq C iqi ( 3 ) 

q e=L,R 

where c\ q (q 9 ) creates (annihilates) an electron with energy e q in lead i = 
L, R. As an initial condition, we employ the grand-canonical ensemble of the 
electrons in the leads at inverse temperature /3 = l/fegT and with electro- 
chemical potentials //l/r- Therefore, the lead electrons are described by the 
equilibrium Fermi function fe(e q ) = {1 + exp[— f3(e q — For the initial 

density matrix, we then have 

£>leads,eq OC exp [-/3(if l eads - fl L N L - fX R N R )] , (4) 

where Ne = J2 q c l q c e q denotes the electron number in the left and right lead, 
respectively. From this follows that all expectation values of the lead operators 
can be traced back to the expression 

( c l'g' c iq) = $U'5 qq 'fe(eq)- (5) 

The terminating sites |1) and \N), the so-called donor and acceptor sites, 
couple via the tunnelling matrix element Ve q to the state \£q) in the respective 
lead. The Hamiltonian describing this interaction has the form 

^wire-lead = ^(^C^Ci + V Rq C Rq C N ) + H.C. (6) 

q 

It will turn out that the influence of the tunnelling matrix elements is com- 
pletely characterised by the spectral density 

r £ (e) = 27r£|Vy 2 5( e - e<? ) (7) 
i 
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which becomes a continuous function of e if the lead modes are dense. If 
all relevant lead states are located in the centre of the conduction band, the 
energy-dependence of the spectral densities is not relevant and can be replaced 
by a constant, I^(e) = T e . This defines the so-called wide-band limit. 



3 Master equation approach 

The computation of stationary currents can be achieved by deriving a master 
equation for the dynamics of the wire electrons. Thereby, the central idea is to 
consider the contact Hamiltonian (6) as a perturbation. From the Liouville-von 
Neumann equation ihg = [H, g] for the total density operator g one obtains 
by standard techniques [29] the approximate equation of motion 

g(t) = - ^[H whc (t) + iJi ea ds, g(t)} 

(8) 

dr [H wire-lead; [-^wire-lead ( ~t) , (>(£)]]■ 

The tilde denotes operators in the interaction picture with respect to the 
molecule and the lead Hamiltonian, X(t) = U^t) X U (t), where U is the 
propagator without the coupling. For the evaluation of Eq. (8) it is essential to 
use an exact expression for the zeroth-order time evolution operator Uo(t). The 
use of any approximation bears the danger of generating artifacts, which, for 
instance, may lead to a violation of fundamental equilibrium properties [29,30]. 

The stationary current defined as the net (incoming minus outgoing) electrical 
current through contact £ is given by minus the time-derivative of the electron 
number in that lead multiplied by the electron charge — e, Ii(t) = e(d/dt){Ne). 
From the master equation (8) follows 

I e (t) = eti[Q(t)N e ] 

g roo , — . . 
= ~ h? JO dr \[^irc-lcad(-7"), [//wire-lead,-^]]), (9) 

where we have used the relation tr A[B,C] = tr[A,B]C. Next, we insert 
the wire-lead Hamiltonian (6), the interaction-picture operator r) = 
Q g exp(ie 9 r) and the expectation values (5). By use of the spectral density 
(7), the remaining sum over the lead states is transformed into an integral 
which in the wide-band limit r^(e) = can be evaluated to read 

h = ^(c[ Cl ) - e^Re fdr J dee^ + ^ h f(e)([ Cl , ~c\(-r)} + ). (10) 

In the following, we specify the master equation (8) and the current formula 
(10) for studying two limiting cases: The first limit U — describes non- 
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interacting electrons. For this situation, we follow the approach of Ref. [31]. 
The second limit is the one of strong Coulomb repulsion in which U is much 
larger than any other energy scale of the problem. Then, only the states with 
at most one excess electron on the wire are relevant. 

In both cases, a diagonal representation of the first term on the right-hand side 
of the master equation (8) is achieved by a decomposition into the eigenbasis 
of the single-particle wire Hamiltonian. In this basis, the fermionic interaction 
picture operators read 

Cn(t) = Y,( n \^)c a e-' lea \ (11) 

a 

where \4> a ) denotes an eigenstate with energy e a . Below, we will need in par- 
ticular the creation and annihilation operators for the sites with direct contact 
to the leads, i.e. \ne) where = 1 and % = N. 



3.1 Non-interacting electrons 



In the limit U = 0, the transport problem defined by the Hamiltonian (1) 
possesses an exact solution which is conveniently derived within a scattering 
approach. However, since one aim of the present work is the comparison of two 
distinct master equations, we only sketch the exact solution for the special case 
of a two- level system in the Appendix A.l and review here the corresponding 
master equation approach [20,31]. 

In general, the relation between the states \<p a ) and the many-particle Hamil- 
tonian (1) is established via the Slater determinant. Alternatively, one can 
resort to Green's functions. In the present case, knowledge of the Green's 
function at time t = is already sufficient. Apart from a prefactor, it is given 
by the expectation value 

Pa P = (4c a ) = P; a . (12) 

Then, one obtains from (10) for the stationary current the expression 

= ^((t>p\ni)(ne\4) a )P a p ~ \(ne\(f>a)\ 2 fe(e a ) 



a 



(13) 



where the index refers to U = 0. It can be shown that the current is inde- 
pendent of the index £, i.e. independent of the contact at which it is evaluated. 
This reflects for a two-probe setting the validity of the continuity equation. 

In order to determine the expectation values P a/ 3, we employ the master equa- 
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tion (8) and obtain for the stationary state the condition 

r r 

i(e a - efs)P al 3 = J2 -w\ {(t>a\ne)(ni\(t)p) [fe(e a ) + fe(ep)] 

-^>2((pa\ne)(n e \(f) a ,) P a > p - ^2{(f>p'\nt)(ne\(j)fi) P a /3> 

a' /3> 

(14) 

In a non-equilibrium situation, the solution of this set of equations generally 
possesses non-vanishing off-diagonal elements, which in some cases turn out 
to be crucial. 



3.2 Strong Coulomb repulsion 



In the limit of strong Coulomb repulsion, U is assumed to be so large that at 
most one excess electron resides on the wire. Thus, the available Hilbert space 
is restricted to the states {|0), c^|0)} Q= i...Ar, which we use for the decomposition 
of the density operator to obtain 

P = |0)poo(0| + £ (4|0)p QO (0| + |0) Poa (0|c a ) + £ 4 |0>p Q/3 (0| C/3 . (15) 

With this ansatz, the current expectation value (10) assumes the form 

Joo = eT e Y^ {^2(4>f3\ne)(n e \(f) a }ffXe a )p a p - \(<j) a \n e }\ 2 f e (e a )p 00 \ , (16) 

a /3 

where / = 1 — /. The decomposition of the master equation (8) into the single- 
particle states cJjO) provides for the stationary state the set of equations 

r r 

e=L,R A 1 

- £ ( <Pa | n t ) (n e | <p a , ) f e (e a , ) p a >p ( 17 ) 

a' 

~ J2i 0/3' I n t ) {nifa ) ft ( ep> ) p a p> 

In order to fully determine the density operator, we need in addition an expres- 
sion for poo which can also be derived from the master equation. A more conve- 
nient alternative is given by the normalisation condition tr p = p 00 + J2 a Paa = 
1. For the sake of completeness, we remark that from the master equation (8) 
follows p a0 = p 0a = in the stationary state. 

It can be shown that if the wire consists of just one site, i.e. for N — 1, both the 
master equation for U = and the one for U = 00 provide identical expressions 
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for the current. The reason for this is that already the Pauli principle inhibits 
the occupation of the molecule by more than one electron. 



3.3 Rotating-wave approximation 



For very weak wire-lead coupling, the coherent time-evolution dominates the 
dynamics of the wire electrons. This means that the largest time-scale of the 
coherent evolution, given by the smallest energy difference, and the dissipative 
time-scale, determined by the coupling rates T^(e), are well-separated, i.e., 

^<|e a -e^| (18) 

for all £ and a ^ (3. Note that this condition is only satisfiable if the energy 
spectrum has no degeneracies. Then for a ^ /3, the master equations (14) 
and (17), which determine the stationary state, are dominated by their left- 
hand side. Consequently, p a p is of the order Y/{e a — ep) such that it can be 
neglected in the limit under consideration. This constitutes the essence of a 
rotating-wave approximation (RWA). The above reasoning is equivalent to the 
assumption that the stationary state is diagonal in the basis of the eigenstates. 
Within such a diagonal ansatz, it is possible to solve both master equations 
analytically and, moreover, to provide a closed expression for the respective 
stationary current. 



3.3.1 RWA for non-interacting electrons 

In the interaction-free case, the stationary state is found by inserting the RWA 
ansatz P a p = P aoi 5 a p into equation (8); after some algebra, we find 

p Wah(e a ) +w*f R (e a ) 

w l + w r ■ ^ 

Thus, the populations are determined by an average over the Fermi functions, 
where the weights 

< = r £ |(n £ |0 Q )| 2 (20) 

are given by the overlap of the eigenstate \<p a ) with the site coupled to lead I. 
Then the average current is readily evaluated to read [32] 

W = eE qn[ [W - fM]. (21) 

This expression represents the limit r — > of the corresponding scattering 
theory 
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3.3.2 RWA for strong Coulomb repulsion 

The corresponding RWA ansatz for the strongly interacting limit reads p a p = 
Paa^a/3- Inserting it into Eq. (17), we find the solution 

J_ Kfde a ) + w%f R (e a ) 
Paa N wtf L (e a )+w*f R (e a y [ZZ) 

with the weight factors w l a defined as above and the normalisation constant 
K r _ , , Wg'hM +w%f R (e a ,) 

The average current follows directly by inserting p aa into (16) and reads 

This current formula differs from the one obtained within RWA for non- 
interacting electrons, c.f. Eq. (21), by the appearance of the normalisation 
factor M and by the Fermi functions / = 1 — / in the denominator. 



4 Open transport channel 



As a first application, we consider a wire for which all on-site energies are at the 
level of the chemical potentials and, moreover, all hopping matrix elements are 
equal, as sketched in Fig. 1. Then, the molecular orbitals are delocalised which 
provides ideal transport along the molecule. For a voltage which is sufficiently 
large such that all molecular orbitals lie within the voltage window, the current 
is in the interaction-free case dominated by the total transmission. Under the 
assumption that all overlaps (20) between molecule eigenstates and sites \n) 
are identical, we find w e a = 1/N. Therefore, the RWA current formula (21) 
becomes 

* = < 25 > 

In particular, we find that the current is independent of the wire length which 
is characteristic for coherent transport of non-interacting electrons through an 
open conduction channel [33]. 

This behaviour is significantly modified by the influence of strong Coulomb 
repulsion. The normalisation factor (23) reads M = N + 1 such that finally 

7 ~ = hW+T)' < 26 > 
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Fig. 2. Stationary current as a function of the wire length for the transport through 
the open channel sketched in Fig. 1 with bias voltage V = 10 A/ h. The other pa- 
rameters are T = 0.1 A and k^T = 0.005A. The dashed lines are a guide to the eye. 




A 



\N) 



MR 



Fig. 3. Level structure of the bridged molecular wire with N = 5 sites. The bridge 
levels are separated by E-q from the donor and acceptor levels |1) and \N). 

For N — 1, this result coincides with Eq. (25), as expected. For a long wire, 
we find that the current decreases oc 1/N, i.e. with the inverse of the wire 
length. This behaviour resembles an ohmic resistor and has been observed in 
the limit of strong Coulomb repulsion [26] also for incoherent hopping of the 
wire electron. 

The numerically computed current beyond RWA is shown in Fig. 2. It fully 
confirms the respective length dependence and, moreover, demonstrates the 
applicability of the rotating-wave approximation in the present case. 



5 Bridged molecular wire 



Next, we consider the bridged molecular wire model sketched in Fig. 3. There, 
the energies of the donor and the acceptor orbitals, |1) and \N), are assumed to 
be close to the chemical potentials of the attached leads, fi L < Ei = E N < /i R . 
The bridge levels E n , n — 2, . . . , N — 1, lie E B 3> A above the chemical 
potential. 

Let us first discuss the eigenstates of the molecule which discern into two 
groups: One group of states is located on the bridge. It consists of iV — 2 levels 
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Fig. 4. Current for the bridged molecular wire model, cf. Fig. 3 comparing the 
non- interacting case with strong Coulomb repulsion. The bias voltage is V = 5 A/ft, 
T = 0.1A, and k B T = 0.005A. The dashed lines mark the analytical results for 
17 = oo. 

with energies in the range [E B — 2 A, E B + 2 A]. In the absence of the driving 
field, these bridge states mediate the super-exchange between the donor and 
the acceptor. The other group consists of the two remaining states. They 
form a doublet whose states are approximately given by (|1) ± |A r ))/y / 2. Its 
splitting can be estimated in a perturbational approach and is approximately 
given by 2A(A/E B ) N ~ 2 [34]. Thus, the wire can be reduced to a two- level 
system with the effective tunnel matrix element A D a = Aexp[— k(N — 2)], 
where k = \n(E B /A). 

The explicit calculation for the two-level system is given in Appendix A. Thus, 
in order to obtain the current for the present case, we just have to replace 
in equations (A. 4) and (A. 13) the tunnel matrix element A by Ada- For 
non-interacting electrons, we find from (A. 4) to lowest order in Ada/T the 
expression 

h = 2 A^l e -MN-2)^ (2?) 

while in the case of strong Coulomb repulsion we employ (A. 13) to obtain 

/oo = e" 2 ^- 2 ). (28) 

In particular, one finds in both cases an exponentially decaying length depen- 
dence of the current [7,15]. Quite remarkably, the strong Coulomb repulsion 
enlarges the current by a factor 2. 

In order to test the quality of the two-level approximation above, we compare 
the analytical result (28) against the numerical solution of the respective mas- 
ter equation. Figure 4 demonstrates the almost perfect agreement between the 
numerical and the analytical solution. Moreover, it conforms the exponentially 
decaying length dependence and the fact that the enhancement of the current 
by a factor 2 owing to Coulomb interaction. For N = 3, the limit Ada T is 
not yet reached which explains the small deviation from expression (28). 
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We close this section with the remark that for bridged molecular wires, the 
rotating-wave approximation derived in Section 3.3 results in I = eT /2h and 
loo = eT/3h. Thus, the RWA even fails to predict qualitatively the observed 
length dependence. The reason for this is that A DA <C T and, thus, the con- 
dition (18) for the applicability of RWA is violated. 



6 Conclusions 



We have derived a master equation approach for the electron transport through 
tight-binding systems in the presence of strong Coulomb repulsion. In contrast 
to prior work, we treat the master equation beyond a rotating-wave approxi- 
mation which extends the range of validity of our approach to intermediately 
strong wire-lead coupling. In particular, bridged molecular wires constitutes 
an example for which our approach provides reliable results while within a 
rotating-wave approximation one obtains qualitatively wrong results. 

With this formalism, we have studied transport properties of two models for 
molecular wires. Thereby, we have worked out the differences to the case of 
non-interacting electrons. A model for which all on-site energies are identi- 
cal, represents a tight-binding version of an open conduction channel. There, 
we find a significant influence of Coulomb repulsion: While in the absence of 
interaction, the current is length independent, it decreases due to Coulomb re- 
pulsion proportional to the wire length. Thus it resembles an ohmic conductor 
even though the transport is fully coherent. 

For the bridged molecular wire model, only the first and the last site have 
energies close to the chemical potentials of the leads, whereas all the other 
sites merely mediate co-tunnelling. We have demonstrated that then the wire 
exhibits the behaviour of a two-level system with an effective tunnel matrix 
element. In particular, we found the surprising effect that Coulomb blocking 
enhances the current by a factor two. 

Comparing the results for the open-channel model and the bridge model, we 
can conclude that the influence of Coulomb repulsion depends sensitively on 
the level structure of the molecule: If many unoccupied molecular orbitals 
have energies close to the chemical potentials of the leads, electron-electron 
interaction reduces the current considerably. 
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A Two-level system 



The bridged molecular wire discussed in Section 5 can be described by a 
conductor that consists of only a donor state |1) and an acceptor state |2), i.e. 
N = 2. In this appendix, we derive explicit results for the transport through 
an unbiased two-level system {E\ — E2 — 0) for Tl = Tr = T and chemical 
potentials such that effectively /l = and /r = 1. These two sites are coupled 
by a tunnel matrix element A. Diagonalising the wire Hamiltonian (2) for 
vanishing Coulomb interaction (U = 0), we obtain the bonding and anti- 
bonding eigenstates and eigenenergies 



-L(|l) + |2», e+ = -A, 
i=(|l>-|2», e_ = A. 



(A.l) 



For weak coupling between the donor and the acceptor, A < T, the master 
equation approach, albeit perturbational in T, still provides the correct be- 
haviour owing to the proper inclusion of off-diagonal elements of the density 
matrix. Within rotating-wave approximation this is no longer the case. 



A . 1 Landauer form 



According to Landauer [35], the coherent transport for non-interacting elec- 
trons can be interpreted as a quantum mechanical scattering process. Thereby 
the in- and outgoing electronic states scattered in the mesoscopic conductor 
are considered as plane waves. As a consequence, the pivotal quantity which 
determines the system's conductance is the transmission probability T(E) and 
the corresponding current can be written in the form 

JdE[f R (E)-f L (E)}T(E). (A.2) 



2nh 
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The transmission can now be calculated via the relation T(E) = T-^TnlGuiE)^ , 
where G(E) — (E — H wirc — iT/2) -1 denotes the retarded Green's function. 
For an unbiased two-level system, we obtain 

T < E) = | (g -.r/2)-Ay < A ' 3) 

Inserting this expression into Eq. (A. 2), one arrives at 



_ er A 2 

~ 2hA 2 + (r/2) 2 ' 

A more explicit calculation can be found, e.g., in Ref. [36]. 



(A.4) 



A. 2 Non-interacting electrons 



For the eigenstates and eigenenergies (A.l), the current formula (13), valid for 
U — 0, reads I Q = (eT/2h) J2 a ,p Pap while the set of equations (14) becomes 

= ^(1-2P ++ ), (A.5) 
= ^(1-2P__), (A.6) 
-2iAP + _ = ^(l-2P + _). (A.7) 

This corresponds to P++ = P__ = 1/2 and P+_ = T/(2r - 4iA) = P*_ + . 
Inserting this solution into equation (13), we obtain for the stationary current 
the result (A.4). 

The quality of the present master equation approach is underlined by the 
fact that it here indeed reproduces even for T ^> A the exact solution. We 
emphasise that this is not the case for the RWA solution (21): Since for this 
approximation by definition P + _ = 0, one obtains the result Io.rwa — eT/2h 
which is independent of the inter-site coupling A. 



A. 3 Strong Coulomb repulsion 



Using the eigenstates and eigenenergies (A.l), one finds that the current with 
Coulomb blocking, Eq. (16), reads 1^ = (eT /2h)J2 a ,/3 Pa/3- Formally, this is 
identical to the corresponding expression for the non-interacting case but with 
the P a p replaced by the matrix elements of the density operator in the basis of 
the single particle states c+|0). The stationary value of these matrix elements 
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is determined by Eq. (14) which for a two- level system reads 



= ^(1 - 2 P++ - p__ - Rep+_), (A.8) 
= ^(1 - 2p__ - P++ - Rep + _), (A.9) 
-2iAp + _ = £(1 - \ P++ - ip__ + p + _), (A.10) 



which corresponds to 



8A 2 + r 2 . . 

P++ = P - = 2(12A 2 + P)' (A - U) 

4Air + r 2 

p+- = p- + = - 2(12A2 + r2) - ( A - 12 ) 



Inserting into Eq. (16), we finally obtain the current 



/ eT (A 13) 

°°~ 2/i3A 2 + (r/2) 2 - 1 ' 
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